Christopher M. Moore
17 August 2016
Cleveland R Users Group
library(package = "deSolve") # Solver for differential equations
library(package = "rootSolve") # Finds roots, Jacobian matrices, and estimates steady-state conditions
library(package = "phaseR") # Tools for graphing
library(package = "MASS") # For 2-D density estimator
library(package = "QPot") # Our package!
library(package = "viridis") # Great color map
The study of the spatial and temporal patterns of the distributions and abundances of organisms, including causes and consequences
Ecology is hierarchically organized:
Today's focus, population ecology, is strongly quantitative and mathematical
\[ \frac{\text{change in a variable}}{\text{change in another variable}} \]
if \( ~x \) is a variable that is dependent on variable \( ~y \), then we can write
\[ \frac{d x}{d y} \]
with time, \( ~t \), as the independent variable
\[ \frac{d x}{d t} \]
Growth of \( ~x \) at a rate, \( ~\alpha \) \[ \frac{dx}{dt} = \alpha x \]
Logistic equation, with growth of \( ~x \) at a rate \( ~\alpha \), and negative density-dependent growth (crowding effect) of strength \( \gamma \) and interaction between \( x \times x \):
\[ \frac{dx}{dt} = \alpha x - \gamma x^2 \]
Competition between species \( ~N_1 \) and \( ~N_2 \)
\[ \frac{dN_1}{dt} = r_1N_1 - \alpha_{11} N_1N_1 - \beta_{21} N_2N_1 \\ \frac{dN_2}{dt} = r_2N_2 - \alpha_{22} N_2N_2 - \beta_{12} N_1N_2 \]
Disease model of susceptible, \( ~S \), infected \( ~I \), and recovered, \( ~R \)
\[ \frac{dS}{dt} = \alpha S - \beta SI + \delta R - d_SS \\ \frac{dI}{dt} = \beta SI -\gamma I - d_II \\ \frac{dR}{dt} = \gamma I -\delta R - d_RR \]
deSolve has set the standard, with:
y0, for each of the variablesThe mathematics of uncontrolled growth are frigtening. A single cell of the bacterium E. coli would, under ideal circumstances, divide every twenty minutes. That is not particularly disturbing until you think about it, but the fact is that bacteria multiply geometrically: one becomes two, two become four, four become eight, and so on. In this way, it can be shown that in a single day, one cell of E. coli could producece a super-colony equal in size and weight to the entire planet earth.—M. Chrichton (1969), The Andromeda Strain (Dell, New York, p. 247)
\[ \frac{dx}{dt} = \alpha x = 2x \]
\[ \frac{dx}{dt} = \alpha x = 2x \]
growth <- function(t, y0, parameters){
with(as.list(c(y0, parameters)), {
dx = alpha*x
list(c(dx))
})
}
growth(y0 = c(x = 1), parameters = c(alpha = 2))
[[1]]
[1] 2
\[ \frac{dx}{dt} = \alpha x = 2x \]
start <- c(x = 1)
parameters <- c(alpha = 2)
hours <- seq(from = 0, to = 24, by = (1/3))
out <- ode(y = start, times = hours, func = growth, parms = parameters)
out[1:3,]
time x
[1,] 0.0000000 1.000000
[2,] 0.3333333 1.947737
[3,] 0.6666667 3.793677
\[ \frac{dx}{dt} = \alpha x = 2x \]
time x
2.400000e+01 7.017045e+20
plot(out[,1], out[,2], xlab = "Time", ylab = "x", type = "l", lwd = 2)
Competition between species \( ~N_1 \) and \( ~N_2 \)
\[ \frac{dN_1}{dt} = r_1N_1 - \alpha_{11} N_1N_1 - \beta_{21} N_2N_1 \\ \frac{dN_2}{dt} = r_2N_2 - \alpha_{22} N_2N_2 - \beta_{12} N_1N_2 \]
comp <- function(t, y0, parameters){
with(as.list(c(y0, parameters)), {
dx <- r1*x - a11*x*x - b21*y*x
dy <- r2*y - a22*y*y - b12*x*y
list(c(dx, dy))
})
}
state <- c(x = 1, y = 4)
parameters <- c(r1 = 1, r2 = 1, a11 = 0.1, a22 = 0.1, b12 = 0.05, b21 = 0.05)
time <- seq(from = 0, to = 25, length.out = 100)
out <- ode(y = state, times = time, func = comp, parms = parameters)
plot(out[,1], out[,2], type = "l", col = "red", xlab = "Time", ylab = "Density of\nspecies 1 (red) and species 2 (blue)", ylim = c(0, max(out[,c(2,3)], na.rm = T)), lwd = 2)
lines(out[,1], out[,3], col = "blue", lwd = 2)
plot(out[,2], out[,3], xlab = "Density of species 1", ylab = "Density of species 2", type = "l", lwd = 2, xlim = c(0, 10), ylim = c(0, 10))
points(out[1, 2], out[1, 3], pch = 16)
plot(out[,2], out[,3], xlab = "Density of species 1", ylab = "Density of species 2", type = "l", lwd = 3, xlim = c(0, 10), ylim = c(0, 10))
for (i in 1:25){
state <- c(x = runif(n = 1, min = 0, max = 10), y = runif(n = 1, min = 0, max = 10))
out <- ode(y = state, times = time, func = comp, parms = parameters)
lines(out[,2], out[,3], lwd = 2)
points(out[1, 2], out[1, 3], pch = 16)
}
comp <- function(t, y0, parameters){
with(as.list(c(y0, parameters)), {
x <- y0[1]
y <- y0[2]
dx <- r1*x - a11*x*x - b21*y*x
dy <- r2*y - a22*y*y - b12*x*y
list(c(dx, dy))
})
}
ff <- flowField(deriv = comp, x.lim = c(0, 10), y.lim = c(0, 10), parameters = parameters, points = 20, add = F, col = "black", arrow.type = "equal")
ff <- flowField(deriv = comp, x.lim = c(0, 10), y.lim = c(0, 10), parameters = parameters, points = 20, add = F, col = "grey50")
for (i in 1:25){
xy.start <- c(runif(n = 1, min = 0, max = 10), runif(n = 1, min = 0, max = 10))
traj <- trajectory(deriv = comp, y0 = xy.start, parameters = parameters, col = "black", t.end = 10, pch = 16, lwd = 2)
}
parameters.ce <- c(r1 = 1, r2 = 1, a11 = 0.1, a22 = 0.1, b21 = 0.25, b12 = 0.25)
ff <- flowField(deriv = comp, x.lim = c(0, 10), y.lim = c(0, 10), parameters = parameters.ce, points = 20, add = F, col = "grey50")
abline(h = 0, v = 0)
for (i in 1:15){
xy.start <- c(runif(n = 1, min = 0, max = 10), runif(n = 1, min = 0, max = 10))
traj <- trajectory(deriv = comp, y0 = xy.start, parameters = parameters.ce, col = "black", t.end = 25, lwd = 2, pch = 16)
}
\[ \frac{dR}{dt} = bR - cRC \\ \frac{dC}{dt} = dRC - eC \]
ff <- flowField(deriv = prey.predator, x.lim = c(0, 3), y.lim = c(0, 3), parameters = parameters, points = 20, add = F, col = "black", xlab = "R", ylab = "C")
\[ \frac{dN_1}{dt} = r_1N_1 - \alpha _{11}N_1^2 + \beta _{21}N_2N_1 \\ \frac{dN_2}{dt} = r_2N_2 - \alpha _{22}N_2^2 + \beta _{12}N_1N_2 \]
ff <- flowField(deriv = mut, x.lim = c(0, 2), y.lim = c(0, 2), parameters = parameters, points = 20, add = F, col = "black")
nc <- nullclines(deriv = mut, x.lim = c(0, 2), y.lim = c(0, 2), parameters = parameters, lwd = 2)
\[ \frac{dN_1}{dt} = r_1N_1 - \alpha _{11}N_1^2 + \beta _{21}N_2N_1 \\ \frac{dN_2}{dt} = r_2N_2 - \alpha _{22}N_2^2 + \beta _{12}N_1N_2 \]
ff <- flowField(deriv = mut, x.lim = c(0, 2), y.lim = c(0, 2), parameters = parameters, points = 20, add = F, col = "black")
nc <- nullclines(deriv = mut, x.lim = c(0, 2), y.lim = c(0, 2), parameters = parameters, lwd = 2)
\[ dN_1 = \left(r_1N_1 - \alpha_{11} N_1N_1 - \beta_{21} N_2N_1\right)dt + \sigma_1 dW_1\\ dN_2 = \left(r_2N_2 - \alpha_{22} N_2N_2 - \beta_{12} N_1N_2\right)dt + \sigma_2 dW_2 \]
parameters <- c(r1 = 1, r2 = 1, a11 = 0.1, a22 = 0.1, b12 = 0.05, b21 = 0.05)
model.deltat <- 0.1
ms <- Model2String(model = comp, parms = parameters, deSolve.form = T, supress.print = T)
ts <- TSTraj(y0 = c(x = 1, y = 1), time = 100, deltat = model.deltat, x.rhs = ms[1], y.rhs = ms[2], sigma = 1, lower.bound = 0)
TSPlot(mat = ts, deltat = model.deltat)
TSPlot(ts, deltat = model.deltat, dim = 2)
TSDensity(mat = ts, deltat = model.deltat, dim = 2, col2d = viridis(100, option = "A"))
The potential relates to the work needed to move from one point to another, with the areas of least potential at surface minima.
The quasi-potential is a tool that yields information about properties of stochastic systems, such as the expected time to escape a basin of attraction, the expected frequency of transitions between basins, and the stationary probability distribution.
Calucualtes the minimal cost path on a surface
Sources for more information
\[ \frac{dN_1}{dt} = r_1N_1 - \alpha_{11} N_1N_1 - \beta_{21} N_2N_1 \\ 0 = r_1N_1 - \alpha_{11} N_1N_1 - \beta_{21} N_2N_1 \\ 0 = N_1\left(r_1 - \alpha_{11} N_1 - \beta_{21} N_2\right) \\ 0 = r_1 - \alpha_{11} N_1 - \beta_{21} N_2 \\ \alpha_{11} N_1 = r_1 - \beta_{21} N_2 \\ N_1^* = \frac{r_1 - \beta_{21} N_2}{\alpha_{11}} \]
ce <- Model2String(model = comp, parms = parameters.ce, deSolve.form = T)
local.1 <- QPotential(x.rhs = ce[1], x.start = 0, x.bound = c(-5, 15), x.num.steps = 2000, y.rhs = ce[2], y.start = 10, y.bound = c(-5, 15), y.num.steps = 2000)
QPContour(surface = local.1, dens = c(500, 500), x.bound = c(-1, 15), y.bound = c(-1, 15))
saddle <- steady(y = c(3,3), func = comp, parms = parameters.ce)
saddle$y
[1] 2.857143 2.857143
global.qp <- QPGlobal(local.surfaces = list(local.1, local.2), unstable.eq.x = c(saddle$y[1]), unstable.eq.y = c(saddle$y[2]), x.bound = c(-5, 15), y.bound = c(-5, 15))
QPContour(surface = global.qp, dens = c(1000, 1000), x.bound = c(-5, 15), y.bound = c(-5, 15), xlim = c(-1, 11), ylim = c(-1, 11))
\[ \frac{dx(t)}{dt} = \alpha x(t)\left(1 - \frac{x(t)}{\beta}\right) - \frac{\delta x^2(t)y(t)}{\kappa + x^2(t)} \\ \frac{dy(t)}{dt} = \frac{\gamma x^2(t)y(t)}{\kappa + x^2(t)} - \mu y^2(t) \]